package scu.maqiang.fes;

import static scu.maqiang.numeric.Constants.GP3;
import static scu.maqiang.numeric.Constants.GW3;
import static scu.maqiang.numeric.MVO.matmul;
import static scu.maqiang.numeric.MVO.mul;
import static scu.maqiang.numeric.MVO.transpose;

import java.util.Arrays;

import scu.maqiang.mesh.MatrixFunc;
import scu.maqiang.mesh.Mesh3H20;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class FES3H202 extends FES {
	public FES3H202(Mesh3H20 mesh) {
		super(mesh, 3);
		eleStiff = FES3H202::elementElasticityStiff;
		eleMass = FES3H202::elementElasticityMass;
		eleSource = FES3H202::elementVolumeForce;
		eleFlux = FES3H202::elementSurfaceForce;
	}

	public static void elementElasticityStiff(int[] ele, double[][] xyz, double[][] ENu, BVPType tp, double[][] K) {
		double[] x = MVO.col(xyz, 0);
		double[] y = MVO.col(xyz, 1);
		double[] z = MVO.col(xyz, 2);
		double[] dx = new double[20];
		double[] dy = new double[20];
		double[] dz = new double[20];
		double[] v = new double[20];
		double xi, et, zt, detJ, wxarea;
		double[] E = ENu[0];
		double[] Nu = ENu[1];
		double mE, mNu;
		double[][] KK = new double[60][60];
		double[][] B = new double[6][60];
		;
		MVO.fill(K, 0.0);
		switch (tp) {
		case COMMON:
			B = new double[6][60];
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					for (int k = 0; k < 3; k++) {
						zt = GP3[k];
						detJ = CG.ISOMap3D(CG::H20ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						mE = MVO.dot_product(v, E);
						mNu = MVO.dot_product(v, Nu);
						double[][] D = Materials.ElasticConstituteMatrix(mE, mNu, tp);
						wxarea = GW3[i] * GW3[j] * GW3[k] * detJ;
						for (int ii = 0; ii < 20; ii++) {// 11, 22, 33, 12, 23, 13
							B[0][3 * ii] = B[3][3 * ii + 1] = B[5][3 * ii + 2] = dx[ii];
							B[1][3 * ii + 1] = B[3][3 * ii] = B[4][3 * ii + 2] = dy[ii];
							B[2][3 * ii + 2] = B[4][3 * ii + 1] = B[5][3 * ii] = dz[ii];
						}
						mul(matmul(transpose(B), matmul(D, B)), wxarea, KK);
						for (int ii = 0; ii < 60; ii++) {
							for (int jj = 0; jj < 60; jj++) {
								K[ii][jj] += KK[ii][jj];
							}
						}
					}
				}
			}
			break;
		case CONSTITUTE_MATRIX_COMMON:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					for (int k = 0; k < 3; k++) {
						zt = GP3[k];
						detJ = CG.ISOMap3D(CG::H20ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						wxarea = GW3[i] * GW3[j] * GW3[k] * detJ;
						for (int ii = 0; ii < 20; ii++) {// 11, 22, 33, 12, 23, 13
							B[0][3 * ii] = B[3][3 * ii + 1] = B[5][3 * ii + 2] = dx[ii];
							B[1][3 * ii + 1] = B[3][3 * ii] = B[4][3 * ii + 2] = dy[ii];
							B[2][3 * ii + 2] = B[4][3 * ii + 1] = B[5][3 * ii] = dz[ii];
						}
						mul(matmul(transpose(B), matmul(ENu, B)), wxarea, KK);
						for (int ii = 0; ii < 60; ii++) {
							for (int jj = 0; jj < 60; jj++) {
								K[ii][jj] += KK[ii][jj];
							}
						}
					}
				}
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}

	public static void elementElasticityMass(int[] ele, double[][] xyz, double[][] r, BVPType tp, double[][] K) {
		double[] x = MVO.col(xyz, 0);
		double[] y = MVO.col(xyz, 1);
		double[] z = MVO.col(xyz, 2);
		double[] dx = new double[20];
		double[] dy = new double[20];
		double[] dz = new double[20];
		double[] v = new double[20];
		double xi, et, zt, detJ, wxarea, c;
		MVO.fill(K, 0.0);
		switch (tp) {
		case CONSIST_MASS:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					for (int k = 0; k < 3; k++) {
						zt = GP3[k];
						detJ = CG.ISOMap3D(CG::H20ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						c = MVO.dot_product(r[0], v);
						wxarea = GW3[i] * GW3[j] * GW3[k] * detJ * c;
						for (int si = 0; si < 20; si++) {
							for (int sj = 0; sj < 20; sj++) {
								K[3 * si][3 * sj] += v[si] * v[sj] * wxarea;
								K[3 * si + 1][3 * sj + 1] += v[si] * v[sj] * wxarea;
								K[3 * si + 2][3 * sj + 2] += v[si] * v[sj] * wxarea;
							}
						}
					}
				}
			}
			break;
		case LUMP_MASS:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					for (int k = 0; k < 3; k++) {
						zt = GP3[k];
						detJ = CG.ISOMap3D(CG::H20ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						c = MVO.dot_product(r[0], v);
						wxarea = GW3[i] * GW3[j] * GW3[k] * detJ * c;
						for (int si = 0; si < 20; si++) {
							K[3 * si][3 * si] += v[si] * wxarea;
							K[3 * si + 1][3 * si + 1] += v[si] * wxarea;
							K[3 * si + 2][3 * si + 2] += v[si] * wxarea;
						}
					}
				}
			}
			break;
		default:
			throw new IllegalArgumentException("Illegal BVPType!");
		}
	}

	public static void elementVolumeForce(int[] ele, double[][] xyz, double[][] f, BVPType tp, double[] eleF) {
		double[] x = MVO.col(xyz, 0);
		double[] y = MVO.col(xyz, 1);
		double[] z = MVO.col(xyz, 2);
		double[] dx = new double[20];
		double[] dy = new double[20];
		double[] dz = new double[20];
		double[] v = new double[20];
		double xi, et, zt, detJ, wxarea, fu, fv, fw;
		Arrays.fill(eleF, 0.0);
		switch (tp) {
		case COMMON:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					for (int k = 0; k < 3; k++) {
						zt = GP3[k];
						detJ = CG.ISOMap3D(CG::H20ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						fu = MVO.dot_product(v, f[0]);
						fv = MVO.dot_product(v, f[1]);
						fw = MVO.dot_product(v, f[2]);
						wxarea = GW3[i] * GW3[j] * GW3[k] * detJ;
						for (int si = 0; si < 8; si++) {
							eleF[3 * si] += v[si] * fu * wxarea;
							eleF[3 * si + 1] += v[si] * fv * wxarea;
							eleF[3 * si + 2] += v[si] * fw * wxarea;
						}
					}
				}

			}
			break;
		default:
			throw new IllegalArgumentException("Illegal BVPType!");
		}
	}

	public static void elementSurfaceForce(int[] ele, double[][] xyz, double[][] f, BVPType tp, double[] eleF) {
		double[] x = MVO.col(xyz, 0);
		double[] y = MVO.col(xyz, 1);
		double[] z = MVO.col(xyz, 2);
		double[] v = new double[8];
		double xi, et, gu, gv, gw, detJ, wxarea;
		Arrays.fill(eleF, 0.0);
		for (int i = 0; i < 3; i++) {
			xi = GP3[i];
			for (int j = 0; j < 3; j++) {
				et = GP3[j];
				detJ = CG.ISOMap3DTo2D(CG::Q4ShapeFunction, x, y, z, xi, et, v);
				gu = MVO.dot_product(f[0], v);
				gv = MVO.dot_product(f[1], v);
				gw = MVO.dot_product(f[2], v);
				wxarea = GW3[i] * GW3[j] * detJ;
				for (int ii = 0; ii < 8; ii++) {
					eleF[3 * ii] += v[ii] * gu * wxarea;
					eleF[3 * ii + 1] += v[ii] * gv * wxarea;
					eleF[3 * ii + 2] += v[ii] * gw * wxarea;
				}
			}
		}
	}

	public void assembleHeatLoad(ScalarFunc[] coef, double[] param, double[][] vecCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(coef, param, vecCoef, FES3H202::elementHeatLoad, tp, RHS);
	}

	public void assembleHeatLoad(double[] coef, double[][] vecCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(coef, vecCoef, FES3H202::elementHeatLoad, tp, RHS);
	}

	public static void elementHeatLoad(int[] ele, double[][] coord, double[][] coef, BVPType tp, double[] eleLoad) {
		double[] x = MVO.col(coord, 0);
		double[] y = MVO.col(coord, 1);
		double[] z = MVO.col(coord, 2);
		double[] dx = new double[20];
		double[] dy = new double[20];
		double[] dz = new double[20];
		double[] v = new double[20];
		double xi, et, zt, detJ, wxarea;
		Arrays.fill(eleLoad, 0.0);
		switch (tp) {
		case COMMON:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					for (int k = 0; k < 3; k++) {
						zt = GP3[k];
						detJ = CG.ISOMap3D(CG::H20ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						wxarea = GW3[i] * GW3[j] * GW3[k] * detJ * MVO.dot_product(v, coef[0]) * MVO.dot_product(v, coef[1]);
						for (int ii = 0; ii < 20; ii++) {
							eleLoad[3 * ii] += dx[ii] * wxarea;
							eleLoad[3 * ii + 1] += dy[ii] * wxarea;
							eleLoad[3 * ii + 2] += dz[ii] * wxarea;
						}
					}

				}
			}
			break;
		case CONSTITUTE_MATRIX_COMMON:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					for (int k = 0; k < 3; k++) {
						zt = GP3[k];
						detJ = CG.ISOMap3D(CG::H20ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						wxarea = GW3[i] * GW3[j] * GW3[k] * detJ;
						double b1 = MVO.dot_product(v, coef[0]);
						double b2 = MVO.dot_product(v, coef[1]);
						double b3 = MVO.dot_product(v, coef[2]);
						double T = MVO.dot_product(v, coef[3]);
						for (int ii = 0; ii < 20; ii++) {
							eleLoad[3 * ii] += dx[ii] * wxarea * b1 * T;
							eleLoad[3 * ii + 1] += dy[ii] * wxarea * b2 * T;
							eleLoad[3 * ii + 2] += dz[ii] * wxarea * b3 * T;
						}
					}

				}
			}
			break;
		default:
			throw new IllegalArgumentException("Illegal BVPType!");
		}

	}
}
